
/*bertrand first order conditions for aids demand.  
s is a nbrnds x 1 matrix of pre-merger revenue shares, order is aunt jemima, hungry jack, log cabin, mrs. butterworth, p.l., aunt jemima light,
hungry jack light, log cabin light, mrs. butterworth light.
p a nbrnds x 1 matrix of pre-merger prices
c is a nbrnds x 1 matrix of pre-merger cost
g_reg is a nbrnds1 x nbrnds1 matrix of coefficients on lnp from the aids model, regular segment, element (i,j) is the slope coefficient of log price i in equation j.
g_light is a nbrnds2 x nbrnds2 matrix of coefficients on lnp from the aids model, light segment, element (i,j) is the slope coefficient of log price i in equation j
b_reg a nbrnds1 x 1 matrix of coefficients on ln(X/P), regular segment
b_light a nbrnds2 x 1 matrix of coefficients on ln(X/P), light segment
w a nbrnds x 1 vector of average revenue shares, overall
w_nest a nbrnds x 1 vector of average revenue shares conditional on nest
x is a 2 x 1 vector of expenditures: reg and light
delta is a 2 x 2 matrix of middle level coefficients on price indexes
beta is a 2x1 vector of middle level coefficients on expenditures
delta_top is overall elasticitity of demand*/


mata
function msyraids(x,s,p,c,g_reg,g_light,b_reg,b_light,w,w_nest,exp,delta,beta,delta_top)
{

	s_reg=s[1..5,.]
	s_light=s[6..9,.]
	p_reg=p[1..5,.]
	p_light=p[6..9,.]
	c_reg=c[1..5,.]
	c_light=c[6..9,.]
	w_reg=w[1..5,.]
	w_light=w[6..9,.]
	w_nest_reg=w_nest[1..5,.]
	w_nest_light=w_nest[6..9,.]
	x_reg=exp[1,.]
	x_light=exp[2,.]
	delta_mid_reg_reg=delta[1,1]
	delta_mid_reg_light=delta[1,2]
	delta_mid_light_reg=delta[2,1]
	delta_mid_light_light=delta[2,2]
	beta_reg=beta[1,.]
	beta_light=beta[2,.]
	
		
	aidfoc=J(9,cols(x),0)
    cept_reg=J(rows(c_reg),1,0)
    cept_light=J(rows(c_light),1,0)

	sh_reg=J(rows(c_reg),cols(x),0)
	sh_light=J(rows(c_light),cols(x),0)
	
	e11=J(1,cols(x),0)
	e22=J(1,cols(x),0)
	e33=J(1,cols(x),0)
	e44=J(1,cols(x),0)
	e55=J(1,cols(x),0)
	e66=J(1,cols(x),0)
	e77=J(1,cols(x),0)
	e88=J(1,cols(x),0)
	e99=J(1,cols(x),0)
	
	e34=J(1,cols(x),0)
	e43=J(1,cols(x),0)
	e89=J(1,cols(x),0)
	e98=J(1,cols(x),0)
	
	e16=J(1,cols(x),0)
	e61=J(1,cols(x),0)
	e27=J(1,cols(x),0)
	e72=J(1,cols(x),0)
	e38=J(1,cols(x),0)
	e83=J(1,cols(x),0)
	e49=J(1,cols(x),0)
	e94=J(1,cols(x),0)
	e48=J(1,cols(x),0)
	e84=J(1,cols(x),0)
	e39=J(1,cols(x),0)
	e93=J(1,cols(x),0)
	
	
	
	/*intercepts for bottom level regular aids equations*/
    for(i=1;i<=rows(c_reg);i++){
        cept_reg[i,1]=s_reg[i,1]+b_reg[i,1]*(w_nest_reg'*ln(p_reg))-b_reg[i,1]*ln(x_reg)-ln(p_reg)'*g_reg[i,.]' 
    }
	
	/*intercepts for bottom level light aids equations*/
	for(i=1;i<=rows(c_light);i++){
        cept_light[i,1]=s_light[i,1]+b_light[i,1]*(w_nest_light'*ln(p_light))-b_light[i,1]*ln(x_light)-ln(p_light)'*g_light[i,.]' 
    }
		
	
    for(j=1;j<=cols(x);j++){
		/*bottom level regular brand share equations*/
		for(i=1;i<=rows(c_reg);i++){
            sh_reg[i,j]=cept_reg[i,1]+b_reg[i,1]*log(x_reg)-b_reg[i,1]*(w_nest_reg'*ln(x[1..5,j]))+ln(x[1..5,j])'*g_reg[i,.]'
        }
		/*bottom level light brand share equations*/
		for(i=1;i<=rows(c_light);i++){
            sh_light[i,j]=cept_light[i,1]+b_light[i,1]*log(x_light)-b_light[i,1]*(w_nest_light'*ln(x[6..9,j]))+ln(x[6..9,j])'*g_light[i,.]'
        }
		/*own elasticities*/
		e11[1,j]=-1+(1+delta_mid_reg_reg)*w_nest_reg[1,1]+beta_reg*(1+delta_top)*w_reg[1,1]+g_reg[1,1]/sh_reg[1,j]+(b_reg[1,1]/sh_reg[1,j])*(delta_mid_reg_reg*w_nest_reg[1,1]+beta_reg*(1+delta_top)*w_reg[1,1])
		e22[1,j]=-1+(1+delta_mid_reg_reg)*w_nest_reg[2,1]+beta_reg*(1+delta_top)*w_reg[2,1]+g_reg[2,2]/sh_reg[2,j]+(b_reg[2,1]/sh_reg[2,j])*(delta_mid_reg_reg*w_nest_reg[2,1]+beta_reg*(1+delta_top)*w_reg[2,1])
		e33[1,j]=-1+(1+delta_mid_reg_reg)*w_nest_reg[3,1]+beta_reg*(1+delta_top)*w_reg[3,1]+g_reg[3,3]/sh_reg[3,j]+(b_reg[3,1]/sh_reg[3,j])*(delta_mid_reg_reg*w_nest_reg[3,1]+beta_reg*(1+delta_top)*w_reg[3,1])
		e44[1,j]=-1+(1+delta_mid_reg_reg)*w_nest_reg[4,1]+beta_reg*(1+delta_top)*w_reg[4,1]+g_reg[4,4]/sh_reg[4,j]+(b_reg[4,1]/sh_reg[4,j])*(delta_mid_reg_reg*w_nest_reg[4,1]+beta_reg*(1+delta_top)*w_reg[4,1])
		e55[1,j]=-1+(1+delta_mid_reg_reg)*w_nest_reg[5,1]+beta_reg*(1+delta_top)*w_reg[5,1]+g_reg[5,5]/sh_reg[5,j]+(b_reg[5,1]/sh_reg[5,j])*(delta_mid_reg_reg*w_nest_reg[5,1]+beta_reg*(1+delta_top)*w_reg[5,1])
		e66[1,j]=-1+(1+delta_mid_light_light)*w_nest_light[1,1]+beta_light*(1+delta_top)*w_light[1,1]+g_light[1,1]/sh_light[1,j]+(b_light[1,1]/sh_light[1,j])*(delta_mid_light_light*w_nest_light[1,1]+beta_light*(1+delta_top)*w_light[1,1])
		e77[1,j]=-1+(1+delta_mid_light_light)*w_nest_light[2,1]+beta_light*(1+delta_top)*w_light[2,1]+g_light[2,2]/sh_light[2,j]+(b_light[2,1]/sh_light[2,j])*(delta_mid_light_light*w_nest_light[2,1]+beta_light*(1+delta_top)*w_light[2,1])
		e88[1,j]=-1+(1+delta_mid_light_light)*w_nest_light[3,1]+beta_light*(1+delta_top)*w_light[3,1]+g_light[3,3]/sh_light[3,j]+(b_light[3,1]/sh_light[3,j])*(delta_mid_light_light*w_nest_light[3,1]+beta_light*(1+delta_top)*w_light[3,1])
		e99[1,j]=-1+(1+delta_mid_light_light)*w_nest_light[4,1]+beta_light*(1+delta_top)*w_light[4,1]+g_light[4,4]/sh_light[4,j]+(b_light[4,1]/sh_light[4,j])*(delta_mid_light_light*w_nest_light[4,1]+beta_light*(1+delta_top)*w_light[4,1])
		
		/*within nest elasticities*/
		e34[1,j]=(1+delta_mid_reg_reg)*w_nest_reg[4,1]+beta_reg*(1+delta_top)*w_reg[4,1]+g_reg[3,4]/sh_reg[3,j]+(b_reg[3,1]/sh_reg[3,j])*(delta_mid_reg_reg*w_nest_reg[4,1]+beta_reg*(1+delta_top)*w_reg[4,1])
		e43[1,j]=(1+delta_mid_reg_reg)*w_nest_reg[3,1]+beta_reg*(1+delta_top)*w_reg[3,1]+g_reg[4,3]/sh_reg[4,j]+(b_reg[4,1]/sh_reg[4,j])*(delta_mid_reg_reg*w_nest_reg[3,1]+beta_reg*(1+delta_top)*w_reg[3,1])
		e89[1,j]=(1+delta_mid_light_light)*w_nest_light[4,1]+beta_light*(1+delta_top)*w_light[4,1]+g_light[3,4]/sh_light[3,j]+(b_light[3,1]/sh_light[3,j])*(delta_mid_light_light*w_nest_light[4,1]+beta_light*(1+delta_top)*w_light[4,1])
		e98[1,j]=(1+delta_mid_light_light)*w_nest_light[3,1]+beta_light*(1+delta_top)*w_light[3,1]+g_light[4,3]/sh_light[4,j]+(b_light[4,1]/sh_light[4,j])*(delta_mid_light_light*w_nest_light[3,1]+beta_light*(1+delta_top)*w_light[3,1])

		/*across nest elasticities*/
		e16[1,j]=(1+b_reg[1,1]/sh_reg[1,j])*(beta_reg*(1+delta_top)*w_light[1,1]+delta_mid_reg_light*w_nest_light[1,1])
		e61[1,j]=(1+b_light[1,1]/sh_light[1,j])*(beta_light*(1+delta_top)*w_reg[1,1]+delta_mid_light_reg*w_nest_reg[1,1])
		
		e27[1,j]=(1+b_reg[2,1]/sh_reg[2,j])*(beta_reg*(1+delta_top)*w_light[2,1]+delta_mid_reg_light*w_nest_light[2,1])
		e72[1,j]=(1+b_light[2,1]/sh_light[2,j])*(beta_light*(1+delta_top)*w_reg[2,1]+delta_mid_light_reg*w_nest_reg[2,1])
		
		e38[1,j]=(1+b_reg[3,1]/sh_reg[3,j])*(beta_reg*(1+delta_top)*w_light[3,1]+delta_mid_reg_light*w_nest_light[3,1])
		e83[1,j]=(1+b_light[3,1]/sh_light[3,j])*(beta_light*(1+delta_top)*w_reg[3,1]+delta_mid_light_reg*w_nest_reg[3,1])
		
		e48[1,j]=(1+b_reg[4,1]/sh_reg[4,j])*(beta_reg*(1+delta_top)*w_light[3,1]+delta_mid_reg_light*w_nest_light[3,1])
		e84[1,j]=(1+b_light[3,1]/sh_light[3,j])*(beta_light*(1+delta_top)*w_reg[4,1]+delta_mid_light_reg*w_nest_reg[4,1])
		
		e49[1,j]=(1+b_reg[4,1]/sh_reg[4,j])*(beta_reg*(1+delta_top)*w_light[4,1]+delta_mid_reg_light*w_nest_light[4,1])
		e94[1,j]=(1+b_light[4,1]/sh_light[4,j])*(beta_light*(1+delta_top)*w_reg[4,1]+delta_mid_light_reg*w_nest_reg[4,1])
		
		e39[1,j]=(1+b_reg[3,1]/sh_reg[3,j])*(beta_reg*(1+delta_top)*w_light[4,1]+delta_mid_reg_light*w_nest_light[4,1])
		e93[1,j]=(1+b_light[4,1]/sh_light[4,j])*(beta_light*(1+delta_top)*w_reg[3,1]+delta_mid_light_reg*w_nest_reg[3,1])
		
		/*first order conditions for each of the 9 brands*/

		aidfoc[1,j]=x_reg*sh_reg[1,j]*((x[1,j]-c_reg[1,1])/x[1,j])*e11[1,j]+x_light*sh_light[1,j]*((x[6,j]-c_light[1,1])/x[6,j])*e61[1,j]+x_reg*sh_reg[1,j]
        

		aidfoc[2,j]=x_reg*sh_reg[2,j]*((x[2,j]-c_reg[2,1])/x[2,j])*e22[1,j]+x_light*sh_light[2,j]*((x[7,j]-c_light[2,1])/x[7,j])*e72[1,j]+x_reg*sh_reg[2,j]
        
	
		aidfoc[3,j]=x_reg*sh_reg[3,j]*((x[3,j]-c_reg[3,1])/x[3,j])*e33[1,j]+x_reg*sh_reg[4,j]*((x[4,j]-c_reg[4,1])/x[4,j])*e43[1,j]+x_light*sh_light[3,j]*((x[8,j]-c_light[3,1])/x[8,j])*e83[1,j]+x_light*sh_light[4,j]*((x[9,j]-c_light[4,1])/x[9,j])*e93[1,j]+x_reg*sh_reg[3,j]

        aidfoc[4,j]=x_reg*sh_reg[4,j]*((x[4,j]-c_reg[4,1])/x[4,j])*e44[1,j]+x_reg*sh_reg[3,j]*((x[3,j]-c_reg[3,1])/x[3,j])*e34[1,j]+x_light*sh_light[3,j]*((x[8,j]-c_light[3,1])/x[8,j])*e84[1,j]+x_light*sh_light[4,j]*((x[9,j]-c_light[4,1])/x[9,j])*e94[1,j]+x_reg*sh_reg[4,j]

        aidfoc[5,j]=x_reg*sh_reg[5,j]*((x[5,j]-c_reg[5,1])/x[5,j])*e55[1,j]+x_reg*sh_reg[5,j]
		
		aidfoc[6,j]=x_light*sh_light[1,j]*((x[6,j]-c_light[1,1])/x[6,j])*e66[1,j]+x_reg*sh_reg[1,j]*((x[1,j]-c_reg[1,1])/x[1,j])*e16[1,j]+x_light*sh_light[1,j]

        aidfoc[7,j]=x_light*sh_light[2,j]*((x[7,j]-c_light[2,1])/x[7,j])*e77[1,j]+x_reg*sh_reg[2,j]*((x[2,j]-c_reg[2,1])/x[2,j])*e27[1,j]+x_light*sh_light[2,j]

        aidfoc[8,j]=x_light*sh_light[3,j]*((x[8,j]-c_light[3,1])/x[8,j])*e88[1,j]+x_light*sh_light[4,j]*((x[9,j]-c_light[4,1])/x[9,j])*e98[1,j]+x_reg*sh_reg[3,j]*((x[3,j]-c_reg[3,1])/x[3,j])*e38[1,j]+x_reg*sh_reg[4,j]*((x[4,j]-c_reg[4,1])/x[4,j])*e48[1,j]+x_light*sh_light[3,j]

        aidfoc[9,j]=x_light*sh_light[4,j]*((x[9,j]-c_light[4,1])/x[9,j])*e99[1,j]+x_light*sh_light[3,j]*((x[8,j]-c_light[3,1])/x[8,j])*e89[1,j]+x_reg*sh_reg[3,j]*((x[3,j]-c_reg[3,1])/x[3,j])*e39[1,j]+x_reg*sh_reg[4,j]*((x[4,j]-c_reg[4,1])/x[4,j])*e49[1,j]+x_light*sh_light[4,j]	
        }
return(aidfoc)
}

end
